A novel method for measuring the bending rigidity of model lipid membranes by simulating tethers 
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The tensile force along a cylindrical lipid bilayer tube is proportional to the membrane's bending modulus and 
inversely proportional to the tube radius. We show that this relation, which is experimentally exploited to mea- 
sure bending rigidities, can be applied with even greater ease in computer simulations. Using a coarse-grained 
bilayer model we efficiently obtain bending rigidities that compare very well with complementary measure- 
ments based on an analysis of thermal undulation modes. We furthermore illustrate that no deviations from 
simple quadratic continuum theory occur up to a radius of curvature comparable to the bilayer thickness. 



INTRODUCTION 

Bilayer-forming lipids are the basic structural component 
of biological cell membranes. In these amphiphilic molecules 
a hydrophilic group is connected to one or two hydrophobic 
hydrocarbon chains. When dissolved into water they spon- 
taneously assemble into a variety of structures. In Nature 
lipid bilayers form the outer plasma membrane of cells as 
well as the walls of the different cellular compartments and 
organelles, such as the endoplasmic reticulum, the Golgi ap- 
paratus, and the nucleus. fT]. 

Lipid bilayer membranes display interesting physics on 
many different length- and time-scales. On atomistic length 
scales this includes questions such as: How do lipid tail length 
and its degree of saturation influence the bilayer state, how 
does a specific hydrophilic head group facilitate solubiliza- 
tion, or how can water permeate the hydrophobic region? 
On somewhat larger scales the embedding of trans-membrane 
proteins or bilayer fusion are being studied. And on scales ex- 
ceeding several times the bilayer thickness, one may ask how 
vesicles are formed and what shape they have, which forces go 
along with a particular bilayer geometry, or how the demis- 
ing of a multicomponent membrane can trigger morphology 
changes. These different sets of questions require different 
techniques for their treatment. In the present article we fo- 
cus on the physics happening on the large scale end, /. e., on 
the continuum level that may be employed on length scales 
beyond a few tens of nanometers, when a membrane may be 
viewed as a twodimensional fluid elastic sheet. 

As is typical in any coarse-graining scheme, many details 
pertaining the a physical system on a given scale get con- 
densed into a few effective parameters on a larger level. In- 
deed, on the continuum level what remains of all lipid de- 
tail are three material parameters - two moduli describing the 
softest deformation, which is bending, and one length scale 
describing a spontaneous curvature. The respective Hamilto- 
nian, proposed in the early 70's iii can be written as a 
surface integral over the entire membrane: 
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Here, the extrinsic curvature K = 1 /Ri + 1 /i?2 is the sum 
of the two local principal curvatures, and the Gaussian curva- 
ture Kc = l/i?ii?2 is their product. The inverse length Cq 



indicates any spontaneous curvature which the bilayer might 
have, so the first term quadratically penalizes deviations of 
the local extrinsic curvature from Cq . The two moduli n and 
K belonging to the two quadratic curvature expressions are 
refeiTed to as bending modulus and saddle splay modulus, re- 
spectively. If the membrane has two identical leaflets, Cq — 
by symmetry, a situation which does seldomly hold for biolog- 
ical membranes but very frequently for artificial lipid bilayers 
and vesicles. Furthermore, since the surface integral over the 
Ricci scalar R can be expressed as a boundary integral plus 
a topological term, the second term in Eqn. ([T]) most often 
only contributes a constant and can then be ignored. Under 
these conditions there remains only a single physical parame- 
ter characterizing the membrane, the bending modulus n, and 
it is thus the most important one to determine. 

Bending rigidities have been measured experimentally by 
various techniques ill SSSSEIiil EH HQ, all ul- 
timately based on one of two general approaches: One may ei- 
ther utiUze the dependence of thermal undulations on a mem- 
brane's rigidity, or measure the force needed to actively bend 
it. The traditional realization of the first approach is to mon- 
itor the fluctuations of vesicles as a function of wavelength 
by light microscopy, a method termed "flicker spectroscopy" 
|5, 6, 7|. A related experimental method is based on micro- 
pipette manipulation techniques. There, the flicker spectrum 
is successively suppressed by increasing the pipette pressure, 
and the bending rigidity can then be obtained from the low- 
tension regime of the tension-area curve iSBHilill- The 
second approach is typically implemented by measuring the 
force needed to pul l nanoscale bilayer tubes (tethers) from 
vesicles fisll . Since the formation of a tube in- 

volves the creation of a high curvature, the work to pull a 
tether is basically done against bending energy, hence the 
modulus K can be determined from it. 

Determination of the bending rigidity is of course equally 
important in computer simulation studies of lipid bilayers, and 
the spectrum of available methods is the same. However, by 
far the most common approach in simulations is flicker spec- 
troscopy, both for atomistic simulations [16, 17, 18] as well 
as for various coarse grained methods H Q IZll H. 
23l I24II . Only recently den Otter and Briels have proposed 



a method by which constraining forces are applied to actively 
deform the membrane [25], and Farago and Pincus have pro- 
posed a scheme based on the change in free energy of deform- 
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ing the bilayer fl^. Unfortunately, both active methods in- 
volve significant technical and conceptual sophistication. This 
may explain why the idea is not commonly employed, despite 
the fact that particularly for stiff membranes fluctuation based 
schemes encounter difficulties (in experiments as well as in 
simulations), because the thermally excited amplitudes de- 
crease with bending modulus and become difficult to resolve 
at some point. 

In the present article we propose an alternative simulation 
approach for studying the curvature elasticity of membranes 
by an active deformation. Our setup essentially involves mea- 
suring the force necessary to hold a membrane tether, and it is 
thus conceptually identical to its experimental "counterpart". 
As we will see, complications of earlier active schemes are 
avoided, and the simulations are very easy to perform and an- 
alyze. We apply this method to a recently proposed coarse- 
grained solvent-free simulation model l22ll23il and find results 
that agree very well with data from the analysis of the thermal 
fluctuations. Moreover, the method permits us to check, up to 
which curvatures the quadratic model from Eqn. ([TJ remains 
valid. Our results indicate that curvature radii close to the bi- 
layer thickness can be imposed without noticeable deviations 
from Eqn. ([T]i. While the precise location for the breakdown 
of quadratic theory may well be model dependent, its valid- 
ity up to length scales comparable to bilayer thickness is in 
agreement with experimental findings [15i1 . 

CURVATURE ELASTICITY 
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FIG. 1: Flicker spectrum of a fluctuating membrane, plot- 
ted in such a way that in the limit q ^ the expression ap- 
proaches K, see Eqn. (|4| The dashed fine is a fit of the form 
(fceT/K + ci (qaY^)^^ which helps to find the asymptotic 
value; the inset shows the unsealed spectrum. The fit leads to 
K = 12.5 kBT, with an error estimated to be ±1 fceT. The 
system is the same lipid bilayer model 
used for the tethers, see below. 
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In this section we first briefly review the fluctuation ap- 
proach towards bilayer elasticity and discuss some of its dif- 
ficulties. We then introduce the alternative scheme based on 
holding a membrane tether. 

Flicker spectroscopy 

The energy expression in Eqn. ([T]i requires knowledge of 
the local membrane curvature. For essentially flat membranes, 
which can be described by specifying their height h{x,y) 
above some reference plane ("Monge-parametrization"), this 
curvature is given by 
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where V is the two-dimensional nabla operator on the base 
plane. The approximation in the second step is the lowest 
order term in a small gradient expansion. On this level the 
Hamiltonian ([T]) becomes quadratic and can be diagonalized 
by Fourier transformation. Assuming an L x L membrane 
patch with periodic boundary conditions, and writing h{r) = 
J2q hq e"' one finds 
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where we for completeness also added a surface tension term, 
S times excess area. From the equipartition theorem we then 
see that the mean squared amplitude of each mode, /. e., the 
fluctuation spectrum or the static structure factor, is given by 
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A fit of the fluctuation spectrum measured in the simulation 
to this expression yields bending modulus k and tension E. 
Since for wave vectors smaller than qmin — ^/S/k the fluc- 
tuations are tension dominated, it is best to simulate at zero 
tension in order to avoid unnecessary damping of the most 
relevant modes. In this case the expression 1 / {q"^ {\hq\'^) L"^) 
approaches k in the limit g — > 0, as is illustrated in Fig.[T]for a 
model simulation B22l 12311 described in more detail below. For 
wave vectors approaching g,nax — 2tt/w, where w is the bi- 
layer thickness, discrete lipid fluctuations such as protrusions 
require a more careful analysis [24]. 

There are clear limitations for the calculation of k using 
thermal fluctuations. The first and obvious one is that large 
values of k lead to very small amplitudes. Considering (i) 
that n/ksT is typically of order 10 and (ii) how strongly the 
amplitudes decay with increasing wave vector, one realizes 
that it requires substantial statistics to be able to resolve the 
spectrum. Particularly unpleasant in this context is that the 
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most important low-q modes equilibrate slowest, with a time 
scale diverging like 

Also, it must be appreciated that the curvatures probed this 
way are very weak. Since (in the tensionless state) (K^) = 



ksT/L^K, the average (root-mean-square) ra- 



dius of curvature is given by i? = (if^)^^/^ ~ k/UbT, 
i. e. several times the box length. These curvatures are much 
weaker than the ones which one usually imposes on systems 
if one studies vesicles, budding, tubes, or fusion. This leaves 
the question open how relevant the measured elastic constants 
really are 



It is also these limitations which den Otter and Briels II25I1 
had in mind when they proposed ways for obtaining larger 
curvatures, either by more advanced sampling techniques, 
such as umbrella sampling, or by explicitly creating undu- 
lations with larger amplitudes by suitable constraints. They 
found that for larger amplitudes the membrane seemed to 
stiffen, which might suggest that the simple curvature elastic 
model underlying Eqn. ([U breaks down. However, an alter- 
native explanation would be provided by the neglect of higher 
order terms in the small gradient approximation for the cur- 
vature in Eqn. (|2]i. Indeed, for a mode h{x,y) = a s'm{qx) 
it is easy to verify that the ratio between linear and nonlinear 
prediction of the curvature energy is given by 



^(qa)-'+0{iqar') , (5) 



which diverges linearly with growing amplitude. While this 
is qualitative what den Otter and Briels observe, the magni- 
tude of their stiffening is bigger than what Eqn. (|5]l would 
predict. The observed deviation for large amplitudes appears 
more likely to be a result of a residual tension stemming from 
their simulations being done at constant box volume. 



Stretching tethers 

Here we present a method for the calculation of k based on 
a different approach. The basic idea is to impose a deforma- 
tion of the membrane, specifically by creating a curved cylin- 
drical vesicle, and then measure the force required to hold it 
in this deformed state. In the experiment such tethers are typ- 
ically created by first attaching adhesive beads to a suitably 
fixated giant vesicle (or a cell) and then pulling it away with 
a laser tweezer that permits the measurement of the involved 
force. In the simulation such a tether can simply be stabi- 
lized by "spanning" a cylindrical vesicle through the simula- 
tion box, across the periodic boundary conditions. One thus 
simulates a system which is perfectly cylindrical (/. e., there 
are no end effects), and the axial pulling force is readily ob- 
tained from the component of the stress tensor along the box- 
spanning direction. 

With a vesicle radius R and a box length in the direction 
of the spanned vesicle, see Fig.|2] the curvature energy is 
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The axial force under the constraint of fixed area A = 2ttRLz 
is obtained from = {dE/dLz)A = 2ttk/R, hence the 
bending modulus is given by 
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Since both and R are easily measurable, k can be readily 
determined in the simulation. In fact, it is this point where im- 
plementing the tether method in a simulation shows its biggest 
advantage over its experimental counterpart: In a real exper- 
iment R can not be measured directly, since its typical mag- 
nitude is sub-optical. It is thus usually re-expressed in terms 
of the membrane tension S, leading to i? = ■\/k/2I] and thus 
K = (F2/27r)^/2E, but then the tension needs to be monitored 
independently by other means. Recently, however, Cuvelier 
et al. ifTsIl devised a clever setup involving two tethers which 
avoids such complications. 

Even though the above analysis is standard in the tether lit- 
erature Ill2lll3lll4ll5ll . it is still only approximate. Notice that 
this time we have neglected thermal fluctuations altogether 
The formula dTJl relies entirely on a "ground state" argument. 
This is justifiable in two ways. First, for not too small radii 
of curvature the fluctuation contribution to the force, as es- 
timated for instance by a simple plane-wave ansatz for the 
cylindrical modes, turns out to be very small. And second, 
the two most obvious effects which fluctuations have on the 
two terms in Eqn. ^ that need to be determined, F^ and R, 
are working in opposite directions. While clearly the mean 
axial force {Fz) will increase (for exactly the same reason 
that it takes a force to pull a fluctuating polymer straight), the 
fluctuation-corrected mean radius (i?) of the vesicle will de- 
crease, since the total area is constant and the area needed 
for fluctuations has to come from somewhere. Within a plane 
wave approximation these two effects cancel. A more accu- 
rate investigation is a fair bit more subtle 112711 . 

By performing various simulations of tethers with differ- 
ent curvature radius R, we can thus address the question how 
far the present quadratic theory remains valid. Assuming 
symmetric membranes, the next terms by which the Hamil- 
tonian density in Eqn. ^ needs to be amended are quartic 
ones, and these are K^, K^Kq, Kq, and the gradient term 
{\7aK){'V°'K), where Va is the metric -compatible covariant 
derivative lEsll . Since for cylinders Kq — and |Va/^|=0, 
the only remaining term is K'^. Adding jK^K'^ to the energy 
density and repeating the steps leading to Eqn. (Q, we then 
find 
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where £4 = ^ K^j n is a characteristic length scale associated 
with corrections beyond quadratic order, and one typically as- 
sumes that it is related to bilayer thickness. 
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MESOSCOPIC MEMBRANE SIMULATION 

To illustrate our method, we have performed mesoscopic 
simulations of a coarse grained lipid bilayer model recently 
developed in our group f22]. Roughly, lipids are represented 
by three consecutive beads of diameter a (our unit length), 
with one end bead being hydrophilic and the two tail beads 
hydrophobic. The latter feature, in addition to an excluded 
volume interaction, an attraction with a tuneable depth e (our 
unit of energy) and range ■ The unit of time is r = u ^Jm/e, 
where m is the unit of mass. By properly choosing Wc and e, a 
wide range of self-assembling fluid bilayer phases of different 
bending rigidities is obtained. More details can be found in 
Refs. Ii22n23il. 

Coarse-grained Molecular Dynamics simulations of a lipid 
systems with Wc = 1.6 a and k^T — 1.1 e were performed 
using the ESPResSo package [29]. The geometries studied 
are summarized in Table U All simulations were performed 
under canonical (NVT) conditions, using a Langevin thermo- 
stat |30] with friction constant F = 1.0 r^^ to keep the tem- 
perature constant. Within a rectangular box with dimensions 
Lx = Ly and L^, using periodic boundary conditions in all di- 
rections, a cylindrical membrane spanning the z-direction was 
initially set up with a radius i?setup chosen in such a way that 
the area per lipid in both leaflets corresponded to the one for 
a flat tensionless bilayer f23\. Upon starting the simulation 
-Rsetup relaxed (typically within about 1000 r) to its equilib- 
rium value (i?), which is smaller than i?setup by about 3-5% 
due to the area required for fluctuations. For this to happen it 
was quite advantageous that the flip-flop rate of lipids between 
the two leaflets is big enough to permit efficient relaxation of 
area-difference strains going along with changes of the mean 
radius. For the integration, a time step of At — 0.005 r was 
used for most of the systems, while in some cases we needed a 
smaller time step of At = 0.001 t in order to obtain accurate 
results. 

RESULTS AND DISCUSSION 

Figure 12] shows two typical snapshots of equilibrated cylin- 
drical vesicles from different simulations. Notice that while 
fluctuations are clearly visible, they are fairly weak, /. e., the 
vesicle is to a very good approximation cylindrical. We use 
the midplane between the two monolayers to denote the aver- 
age radius (R). It is determined by first identifying the axis, 
next finding the average distance of the second tail bead of the 
outer leaflet to this axis, i?out, and the equivalent for the inner 
leaflet, Rin. We then take (i?) ^{Rout + Rin), where the 
average is taken during long production runs typically extend- 
ing over 10000-20000 r. Errors are determined via a blocking 
analysis. During these runs we also measure the stress tensor 
aij using the virial theorem |31]. Figure[3]shows a typical ex- 
ample of the running average of the three normal stress com- 
ponents. As we observe, the axial stress cr^z has a finite value. 
In contrast, and ayy (as well as all off-diagonal compo- 




FIG. 2: Snapshots of two tether simulations with 20 000 lipids 
and different radii of curvature: a) (i?) = 24crb) (R) = 12 a. 
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FIG. 3: Running average of the diagonal components axx, 
ayy and a^z of the stress tensor for a cylindrical vesicle with 
(i?) = 70(7. 

nents not shown here) approach zero. This is expected, since 
no stress is being transported across the x- and y-direction. 
The error in azz is also determined via a blocking analysis. 

One more point concerning the calculation of the stress ten- 
sor should be mentioned. In general, deriving accurately the 
stress (or the pressure) from molecular dynamics simulations 
is not a trivial aspect. Stress is a collective property with high 
statistical uncertainty owing to the fluctuations of the instan- 
taneous configurations. In common simulations these fluctu- 
ations are very high due to the relatively small size (number 
of particles) of the systems studied. Therefore large systems 
and/or long simulation runs are needed. This point is partic- 
ularly sever in the present case since, as visible in Figure [3] 
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40 



50 
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90 
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5000 
10000 
20000 



24.0±0.5 
3.96±0.9 



47.8±0.6 
1.8±0.7 



16.1±0.4 
5.4±0.8 



12.2±0.3 
6.12±0.8 



24.0±0.4 
3.24±0.7 



47.7±0.6 
1.7±0.8 



9.8±0.2 
8.28±0.8 



8.3±0.2 
10.44±0.8 



16. ±0.4 
4.68±0.8 



7.2±0.1 
11.52±0.8 



6.3±0.1 
12.6±0.8 



12.2±0.2 
6.8±0.8 



24.0±0.4 
3.5±0.8 



5.7±0.1 
13.68±0.8 



9.9±0.2 
8.28±0.8 



8.3±0.1 
10.5±0.7 



16.0±0.3 
4.7±0.8 



6.3±0.1 
12.5±0.8 



12.1±0.2 9.9±0.2 8.3±0.1 
6.6±0.8 8.3±0.8 10.4±0.8 



TABLE I: List of all simulation geometries, sorted by the number of lipids (first column) and their tether length (in units of a, 
first row). In all tuples the top value indicate (i?) (in units of a), the bottom value is (F^) (in units of e/cr). We always used 
= 1-6 a and k^T — 1.1 e. 
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FIG. 4: Tensile force {Fz) as a function of cylinder radius 
(i?) for the systems with N = 5000 lipids (with the ex- 
ception of the open symbol, which has iV = 10 000 lipids). 
The solid line is a fit to Eqn. (Q, leading to k = 11.7 /cbT. 
The inset shows the combination FzR/2ti as a function of 
rescaled curvature w/{R), and the solid line again indicates 
K ^ 11.7 kBT. 



we need to determine very small values of the stress. We have 
found that, in order to obtain reliable values, the common time 
step used in coarse grained simulations. At = 0.01 r, is too 
long. With this choice (Txx and a-yy approached values which 
were significantly different from zero, a clear sign of a sys- 
tematic integration error. We thus used the smaller time step 
At — 0.005 T, and for the cases of very small forces, /. e. large 
{R), even a time step of At = 0.001 t. 

The tensile force in z direction is obtained from the stress 
via = (TzzLxLy. Figure |4] shows the force for the sys- 
tems with N = 5000 lipids as a function of the average cylin- 



der radius (R). As the radius increases, the tensile force F^ 
decreases in accord with Eqn. O. This is seen even better 
by looking at the combination FzR/2'k, which is shown in 
the inset of Fig.|4]as a function of rescaled curvature w/ {R). 
Notice that within the error bars of our simulation this ex- 
pression is perfectly compatible with a constant; a fit gives 
n/k-QT = 11.7 ± 0.2, in very good agreement with the value 
n/k-QT — 12.5 ± 1 obtained from an analysis of thermal un- 
dulation modes (see Fig. [T]!. A possible quadratic deviation, 
as suggested by Eqn. (O, can not be identified with any statis- 
tical significance. This is all the more amazing when we see 
that the most strongly curved cylinder has w/ R k, 0.9, /. e., a 
radius of curvature which is only 10% larger than the bilayer 
thickness w\ Stated differently, the length from Eqn. (O 
must be a fair bit smaller than w. This remarkable robustness 
of the simple quadratic Helfrich theory down to such small 
radii of curvature might of course be a special feature of the 
particular model we have studied, and it would be worthwhile 
to subject other coarse-grained lipid models to a similar test. 
But the fact that extremely high curvatures can be imposed 
without noticing deviations from quadratic continuum theory 
is in accord with common practice in tether pulling experi- 
ments, where the radii of curvature of these membrane tubes 
are typically in the 10-40 nm range, apparently without ever 
having triggered the need to include higher order corrections 
to the elastic behavior llisll . 

Another practical aspect of our proposed method, is related 
to numerical efficiency. The traditional method of analyzing 
the thermal fluctuation spectrum requires very long simula- 
tion runs, since (i) large systems need to be studied in order 
to have a series of wave vectors in the regime where contin- 
uum methods are applicable, and (ii) these long wavelength 
modes take a particularly long time to equilibrate (the relax- 
ation time of bending modes scales with the fourth power of 
wavelength). Strictly speaking our method also requires large 
systems to be studied in order to extrapolate to the zero cur- 
vature limit (/. e., R oo). However, as we have seen in 
Fig.m the asymptotic limit in our case is already reached for 
fairly small systems which still have a significant curvature. 
If the same holds for other lipid models - a fact that needs 
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FIG. 5: Tensile force Fz as a function of system size for sys- 
tems of different ratio LJN: •: 0.01, O: 0.008, ■: 0.006, □: 
0.004. The corresponding average radius (i?) is als indicated. 
A slight horizontal offset is included to improve visibility of 
the error bars. 

to be checked - their value of k can also be determined via 
the tether method using fairly small systems and correspond- 
ing small simulation times. For example, the point in Figure 
|4]having w/{R) = 0.5 is taken from a run of only 2-3 days 
on a single AMD Opteron 2.2GHz processor. For the same 
system the analysis of the thermal fluctuations needs at least 
one month on the same machine. 

As a final point we would like to address the issue of finite 
size effects. While the tether force obviously depends on its 
radius, Eqn. (|7]l suggests that tether length (or equivalently, 
the number N of lipids) is irrelevant. This is indeed rigor- 
ously true in the ground state, but fluctuations might change 
the picture. We have thus repeated our simulations for sys- 
tems with a different number of lipids and checked, whether 
systems with a fixed ratio Lz/N (i. e., essentially identical ra- 
dius) but varying N show any noticeable systematic change 
in the tensile force Fz- Figure |5] shows the results of such 
simulations. As can be seen, the measured forces are, at least 
within our error bars, compatible with a constant value for a 
fixed ratio Lz/N. No finite size effect is detectable. Notice 
that this also provides another independent check that fluctu- 
ations in our system, even though present, are a subdominant 
effect compared to the main "signal" which is well described 
by ground state theory. 



CONCLUSIONS 

We have presented a new method for calculating the bend- 
ing rigidity of lipid membranes in simulations. It involves the 



simulation of cylindrical membrane tethers, spanned across 
the periodic boundary conditions of the simulation box, and 
measuring their equilibrium radius as well as the tensile force 
they exercise on the box. In contrast to fluctuation based 
schemes, which monitor thermally excited shape deforma- 
tions, our approach actively imposes a deformation on the sys- 
tem and measures the restoring force and is thus not limited 
to the regime of deformations accessible by thermal energy. 
In fact, thermal undulations only contribute a small correc- 
tion to the main observable, in stark contrast to fluctuation 
schemes in which they provide the dominant signal. For this 
reason our method is very efficient, also applicable to stiff 
membranes which show very small undulations to begin with, 
and does not crucially depend on the relaxation of very slow 
long wavelength modes. The straightforward access to strong 
bending permits a check of quadratic continuum theory, with- 
out running into difficulties of Monge gauge and its lineariza- 
tion. For the coarse-grained lipid model we explicitly stud- 
ied we showed continuum theory to be applicable up to cur- 
vatures comparable to bilayer thickness. Finite size effects 
would originate from fluctuations and are thus also weak; in 
our runs they were not detectable. 

We believe that this method provides a powerful alternative 
to the existing schemes that is worth to be applied to other 
existing coarse-grained models. Not only is an independent 
measurement of the elastic modulus very valuable, determin- 
ing the range of validity of continuum theory for each model 
would be an important bit of knowledge, given that the cur- 
vatures that are regularly imposed in simulations exceed ther- 
mally excited ones by at least one or two orders of magnitude. 
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